%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 粒子群算法
% 简介：粒子群算法寻找最大值
% 作者：Zhaojiang
% 日期：2023/10/10
% 企鹅：277746470
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all;clear;clc;
%% target fucntion
[x,y]=meshgrid(-4:0.1:4);
z=target(x,y);
mesh(x,y,z);
title('Target Function')
xlabel('x');ylabel('y');zlabel('z');
pause(2);
close all;clear;clc;
%% 设定算法参数
rng(1234);
var_num = 2;% 2个求解变量
lower_bound = -4;% 解的最小值（小数后4位）
upper_bound = 4;% 解的最大值（小数后4位）
bound_length = upper_bound-lower_bound;
max_iteration = 200;% 最大寻找次数
particle_size = 100;% 粒子群大小
w = 1;% 惯性因子
c1 = 2;% 个体学习因子
c2 = 2;% 社会学习因子
v_max = 1;% 粒子最大速度
%% 算法数据初始化
x = lower_bound+bound_length*rand(particle_size,var_num);% 粒子位置
v = 0.2*v_max*rand(particle_size,var_num);% 粒子速度
personnal_x = x;% 个体位置
personnal_fitness = target(x(:,1),x(:,2));% 个体适应度
[global_fitness,i] = min(personnal_fitness);% 社会位置
global_x = personnal_x(i,:);% 社会适应度
% best_fitness = global_fitness;
best_iteration = 0;
all_best_fitness =zeros(1,max_iteration);% 适应度记录
%% 开始迭代求解
for iteration=1:max_iteration
    % 更新位置信息
    for j = 1:particle_size
        for k = 1:var_num
            v(j,k) = w*v(j,k)+c1*rand*(personnal_x(j,k)-x(j,k))...
                     +c2*rand*(global_x(k)-x(j,k));
            if v(j,k) >=v_max
                v(j,k) = v_max;
            elseif v(j,k) <= -v_max
                v(j,k) = -v_max;
            end
            x(j,k) = x(j,k) + v(j,k);
        end
    end
    % 更新个体和社会信息
    for j = 1:particle_size
        fitness = target(x(j,1),x(j,2));
        if fitness > personnal_fitness(j)
            % 社会信息更新
            if fitness > global_fitness
                global_x = x(j,:);
                global_fitness = fitness;
                best_iteration = iteration;
            end
            % 个体信息更新
            personnal_x(j,:) = x(j,:);
            personnal_fitness(j) = fitness;
        end
    end
    all_best_fitness(iteration) = global_fitness;
end
%% 处理结果
% 打印结果
max_x = global_x(1);
max_y = global_x(2);
max_ans = target(max_x,max_y);
disp 最优迭代次数
disp(best_iteration)
disp 最大值位置
fprintf("x=%f\ny=%f\nz=%f\n",max_x,max_y,max_ans);
% 可视化结果
figure('position',[200 200 1280 480])
subplot(1,2,1);hold on;
[x,y]=meshgrid(-4:0.1:4);
z=target(x,y);
mesh(x,y,z);
title('Target Function')
xlabel('x');ylabel('y');zlabel('z');
view(-37.5,30)
pause(0.1)
scatter3(max_x,max_y,max_ans,'r*');
subplot(1,2,2);hold on;
x = 1:max_iteration;
y1= all_best_fitness;
try
    assert(length(x)==length(y1))
    plot(x,y1,'-','LineWidth',1.2)
    xlabel('x');ylabel('y');
    title('迭代数据')
catch
    disp("plot failed")
    disp(size(x))
    disp(size(y1))
end